#include "linear.h"
#include <cilk.h>

int32_t
run_lcs_linear(char* s1, uint32_t sz1, char* s2, uint32_t sz2, char **lcs)
{
  // End of s1/s2?
  if (sz1 == 0 || sz2 == 0) {
    *lcs = NULL;
    return 0;
  }

  uint32_t i, j;

  // Always make s1 smaller string
  if (sz1 > sz2) {
    uint32_t sztemp = sz2;
    char *stemp = s2;
    s2 = s1;
    sz2 = sz1;
    s1 = stemp;
    sz1 = sztemp;
  }

  // s1 is unit string?
  if (sz1 == 1) {
    for (j = 0; j < sz2; j++) {
      if (s1[0] == s2[j]) {
	*lcs = (char *) malloc(sizeof(char));
	(*lcs)[0] = s1[0];
	return 1;
      }
    }

    *lcs = NULL;
    return 0;
  }
  
  uint32_t k;
  char *lcs_b, *lcs_e;
  int32_t *len_b, *len_e;
  int32_t max, l_b, l_e;

  // Split s1 into two substrings at i-th char
  i = sz1 / 2;

  l_b = cilk_spawn lcs_len(s1, i, s2, sz2, &len_b, 0);
  l_e = lcs_len(s1 + i, sz1 - i, s2, sz2, &len_e, 1); // reverse
  cilk_sync;

  if (l_b < 0 || l_e < 0) {
    return -1;
  }

  max = -1;
  k = 0;
  for (j = 0; j < sz2 + 1; j++) {
    if ((len_b[j] + len_e[sz2 - j]) > max) {
      max = len_b[j] + len_e[sz2 - j];
      k = j;
    }
  }

  free(len_b);
  free(len_e);

  l_b = cilk_spawn run_lcs_linear(s1, i, s2, k, &lcs_b);
  l_e = run_lcs_linear(s1 + i, sz1 - i, s2 + k, sz2 - k, &lcs_e);
  cilk_sync;

  if (l_b < 0 || l_e < 0) {
    return -1;
  }

  *lcs = (char *) malloc(sizeof(char) * (l_b + l_e));
  
  if (!*lcs) {
    fprintf(stderr, "Out of memory!\n");
    return -1;
  }

  memcpy(*lcs, lcs_b, l_b);
  memcpy(*lcs + l_b, lcs_e, l_e);

  free(lcs_b);
  free(lcs_e);
  
  return l_b + l_e;
}

int32_t
lcs_len(char *s1, uint32_t sz1, char *s2, uint32_t sz2, int32_t **len, 
	char rev)
{
  uint32_t i, j, size;
  int32_t *prev, *cur;

  size = sizeof(int32_t) * (sz2 + 1);
  cur = (int32_t *) malloc(size);
  prev = (int32_t *) malloc(size);

  if (!cur || !prev) {
    fprintf(stderr, "Out of memory!\n");
    return -1;
  }

  memset(cur, 0, size); // Zero out *len

  if (!rev) {
    for (i = 0; i < sz1; i++) {
      memcpy(prev, cur, size); // prev = cur[:]
      for (j = 0; j < sz2; j++) {
	if (s1[i] == s2[j]) {
	  cur[j + 1] = prev[j] + 1;
	}
	else {
	  cur[j + 1] = MAX(cur[j], prev[j + 1]);
	}
      }
    }
  }
  else {
    for (i = 0; i < sz1; i++) {
      memcpy(prev, cur, size); // prev = curr[:]
      for (j = 0; j < sz2; j++) {
	if (s1[sz1 - 1 - i] == s2[sz2 - 1 - j]) {
	  cur[j + 1] = prev[j] + 1;
	}
	else {
	  cur[j + 1] = MAX(cur[j], prev[j + 1]);
	}
      }
    }
  }

  free(prev);
  *len = cur;

  return sz2 + 1;
}
